if (!require("pacman"))
install.packages("pacman")
pacman::p_load(tidyverse,
dplyr,
ggplot2,
tidyr,
scales,
readr,
ggthemes,
rnaturalearth,
sf,
leaflet
)Question 2
# Loading the main datasets for visualisations
livestock_data <- read.csv("./data/Production_Crops_Livestock_E_All_Data/Production_Crops_Livestock_E_All_Data.csv")
area_codes <- read.csv("./data/Production_Crops_Livestock_E_All_Data/Production_Crops_Livestock_E_AreaCodes.csv")
elements <- read.csv("./data/Production_Crops_Livestock_E_All_Data/Production_Crops_Livestock_E_Elements.csv")
item_codes <- read.csv("./data/Production_Crops_Livestock_E_All_Data/Production_Crops_Livestock_E_ItemCodes.csv")# Identifying livestock commodities that we want
livestock_items <- c("Beef and Buffalo Meat, primary", "Fat of pigs", "Meat, Poultry", "Sheep and Goat Meat", "Milk, Total")
# Filtering for livestock production and trade data
livestock_data_filtered <- livestock_data %>%
filter(Item %in% livestock_items,
Element %in% c("Production", "Area harvested", "Yield")) %>%
select(Area, Item, Element, starts_with("Y"))
columns_to_pivot <- grep("^Y[0-9]{4}$", names(livestock_data_filtered), value = TRUE)
# Reshaping the data to long format
livestock_data_long <- livestock_data_filtered %>%
pivot_longer(
cols = all_of(columns_to_pivot),
names_to = "Year",
values_to = "Value"
) %>%
mutate(Year = as.numeric(gsub("Y", "", Year))) # To convert Year to numeric
# Remove rows with missing Year values if there are any
livestock_data_long <- livestock_data_long %>%
filter(!is.na(Year))
# Aggregate data by year and element
aggregated_data <- livestock_data_long %>%
group_by(Year, Element) %>%
summarise(Total_Value = sum(Value, na.rm = TRUE), .groups = "drop")
# Step 4: Decode area names and item names
area_decoded <- area_codes %>%
rename(Area_Code = `Area.Code`, Area_Name = Area)
item_decoded <- item_codes %>%
rename(Item_Code = `Item.Code`, Item_Name = `Item`)
# For plot 2
commodity_trends <- livestock_data_long %>%
filter(Element == "Production") %>%
group_by(Year, Item) %>%
summarise(Total_Production = sum(Value, na.rm = TRUE), .groups = "drop")
# For plot 3
# Load Natural Earth data for country boundaries
world <- ne_countries(scale = "medium", returnclass = "sf")
# Prepare data for map-based visualizations
# Summarize livestock production by country
map_data <- livestock_data_long %>%
filter(Element %in% c("Production", "Export Quantity", "Import Quantity")) %>%
group_by(Area, Element) %>%
summarize(Total_Value = sum(Value, na.rm = TRUE)) %>%
ungroup()
# Merge map data with country boundaries
map_data <- map_data %>%
left_join(world, by = c("Area" = "admin"))
# Create separate datasets for production, export, and import
map_production <- map_data %>%
filter(Element == "Production")
map_trade <- map_data %>%
filter(Element %in% c("Export Quantity", "Import Quantity")) %>%
mutate(Trade_Type = ifelse(Element == "Export Quantity", "Export", "Import"))
#For Interactive plot 4
# Ensure `geometry` column is in correct SF format
map_production_sf <- map_production %>%
st_as_sf()
population_data <- data.frame(
Area = unique(map_production$Area),
Population = runif(length(unique(map_production$Area)), 1e6, 1e8) # Simulated population
)
# Merge population data to calculate per capita production
map_production <- map_production %>%
left_join(population_data, by = "Area") %>%
mutate(Per_Capita_Production = Total_Value / Population)
# Convert to sf object
map_production_sf <- map_production %>%
st_as_sf()
# Step 2: Define color palette
production_pal <- colorNumeric("YlGn", domain = map_production_sf$Total_Value)
per_capita_pal <- colorNumeric("Blues", domain = map_production_sf$Per_Capita_Production)# Load required libraries
library(ggplot2)
library(dplyr)
# Filter the data to only include "Production" for the Element
aggregated_data_production <- aggregated_data %>%
filter(Element == "Production")
# Create the line plot for production only
ggplot(data = aggregated_data_production, aes(x = Year, y = Total_Value / 1e6, color = Element)) +
geom_line(linewidth = 1.2) +
labs(
title = "Global Trends in Livestock Production",
x = "Year",
y = "Volume (Millions of Tons)",
color = "Element"
) +
theme(
plot.title = element_text(size = 22, face = "bold", hjust = 0.5),
axis.title.x = element_text(size = 19),
axis.title.y = element_text(size = 19),
axis.text = element_text(size = 12),
plot.caption = element_text(size = 14, hjust = 1, color = "grey30"),
legend.position = "right", # Legend on the right
legend.direction = "vertical", # Stack legend vertically
legend.title = element_text(size = 19),
legend.text = element_text(size = 18),
panel.background = element_rect(fill = "white", color = NA), # Set white background
plot.background = element_rect(fill = "white", color = NA), # Ensure white plot background
panel.grid.major = element_line(color = "grey90"), # Subtle grid lines
panel.grid.minor = element_blank() # Remove minor grid lines
) +
scale_y_continuous(labels = scales::comma)# Plot 2: Trends in individual livestock commodities (Production Only)
ggplot(data = commodity_trends, aes(x = Year, y = Total_Production / 1e6, fill = Item)) +
geom_area(position = "stack", alpha = 0.7) +
labs(
title = "Trends in Production of Livestock Commodities",
x = "Year",
y = "Production Volume ( Millions - Tons )",
fill = "Commodity"
) +
theme(
plot.title = element_text(size = 22, face = "bold", hjust = 0.5),
axis.title.x = element_text(size = 19),
axis.title.y = element_text(size = 19),
axis.text = element_text(size = 12),
plot.caption = element_text(size = 14, hjust = 1, color = "grey30"),
legend.position = "right", # Legend on the right
legend.direction = "vertical", # Stack legend vertically
legend.title = element_text(size = 19),
legend.text = element_text(size = 18),
panel.background = element_rect(fill = "white", color = NA), # Set white background
plot.background = element_rect(fill = "white", color = NA), # Ensure white plot background
panel.grid.major = element_line(color = "grey90"), # Subtle grid lines
panel.grid.minor = element_blank() # Remove minor grid lines
)# Plot 3: Livestock Production on a World Map
ggplot() +
geom_sf(data = world) +
geom_sf(data = map_production, aes(geometry = geometry, fill = Total_Value), color = "black") +
scale_fill_viridis_c(name = "Production (tons)", option = "C") +
scale_fill_gradient(
low = "khaki1",
high = "green4"
)+
labs(
title = "Global Livestock Production",
subtitle = "Production volumes (tons) by country",
caption = "Source: FAO Livestock Data"
) +
theme_minimal()# Enhanced Leaflet Map
leaflet() %>%
# Add a base map
addProviderTiles(providers$CartoDB.Positron) %>%
# Add Total Production Layer
addPolygons(
data = map_production_sf,
fillColor = ~production_pal(Total_Value),
color = "black",
weight = 1,
fillOpacity = ~pmax(Total_Value / max(map_production_sf$Total_Value, na.rm = TRUE), 0.2), # Dynamic opacity
label = ~paste(
"<strong>Country:</strong> ", Area, "<br>",
"<strong>Total Production:</strong> ", scales::comma(Total_Value), " tons"
),
popup = ~paste(
"<strong>Country:</strong> ", Area, "<br>",
"<strong>Total Production (tons):</strong> ", scales::comma(Total_Value), "<br>",
"<strong>Population:</strong> ", scales::comma(Population), "<br>",
"<strong>Per Capita Production (tons):</strong> ", round(Per_Capita_Production, 2)
),
highlightOptions = highlightOptions(
weight = 3,
color = "blue",
bringToFront = TRUE
),
group = "Total Production"
) %>%
# Add Per Capita Production Layer
addPolygons(
data = map_production_sf,
fillColor = ~per_capita_pal(Per_Capita_Production),
color = "black",
weight = 1,
fillOpacity = ~pmax(Per_Capita_Production / max(map_production_sf$Per_Capita_Production, na.rm = TRUE), 0.2), # Dynamic opacity
label = ~paste(
"<strong>Country:</strong> ", Area, "<br>",
"<strong>Per Capita Production:</strong> ", round(Per_Capita_Production, 2), " tons"
),
popup = ~paste(
"<strong>Country:</strong> ", Area, "<br>",
"<strong>Per Capita Production (tons):</strong> ", round(Per_Capita_Production, 2), "<br>",
"<strong>Total Production (tons):</strong> ", scales::comma(Total_Value), "<br>",
"<strong>Population:</strong> ", scales::comma(Population)
),
highlightOptions = highlightOptions(
weight = 3,
color = "green",
bringToFront = TRUE
),
group = "Per Capita Production"
) %>%
# Add Layer Control
addLayersControl(
baseGroups = c("Base Map"),
overlayGroups = c("Total Production", "Per Capita Production"),
options = layersControlOptions(collapsed = FALSE)
) %>%
# Add Legend for Total Production
addLegend(
pal = production_pal,
values = map_production_sf$Total_Value,
title = "Total Production<br>(in tons)",
position = "bottomleft",
group = "Total Production"
) %>%
# Add Legend for Per Capita Production
addLegend(
pal = per_capita_pal,
values = map_production_sf$Per_Capita_Production,
title = "Per Capita<br>Production (tons)",
position = "bottomright",
group = "Per Capita Production"
)